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We bring the EKRT framework, which combines perturbative QCD (pQCD) minijet production 
with gluon saturation and hydrodynamics, to next-to-leading order (NLO) in pQCD as rigorously 
as possible. We chart the model uncertainties, and study the viability and predictive power of the 
model in the light of the RHIC and LHC measurements in central A+A collisions. In particular, 
we introduce a new set of measurement functions to define the infrared- and collinear-safe minijet 
transverse energy, Et, in terms of which we formulate the saturation. We update the framework 
with the EPS09 NLO nuclear parton distributions (nPDFs), and study the propagation of the nPDF 
uncertainties into the computed Et, saturation scales and the final-state multiplicities. The key 
parameters, which need to be fixed using the measurements, are identified, and their correlation is 
discussed. We convert the saturated minijet Et into QCD-matter initial conditions for longitudinally 
boost-invariant ideal hydrodynamics. We compute the charged-particle multiplicities and identified 
bulk hadron pt spectra in 5% most central Au+Au collisions at RHIC and Pb+Pb at the LHC. We 
obtain an encouragingly good agreement with the experimental data, simultaneously at RHIC and 
LHC, showing that the approach has a definite predictive power. 

PACS numbers: 24.85.+p, 12.38.Bx, 25.75.Nq, 24.10.Nz 



I. INTRODUCTION 

Different models based on the phenomenon of gluon 
saturation have offered appealing frameworks to predict 
and explain the observed charged-particle multiplicities 
in Pb+Pb collisions at v / ijvw = 2.76 GeV at the Large 
Hadron Collider (LHC), as well as their systematics from 
the Relativistic Heavy Ion Collider (RHIC) to LHC: for 
a collection of the LHC predictions, see Ref. jl|, a com- 
pilation in Ref. [2[ , and a comparison with the first AL- 
ICE data in Ref. 0. Among these is also the so-called 
"EKRT" model introduced in Ref. which combines 
perturbative QCD minijet production |5[ with a satura- 
tion of gluons 1 to compute (instead of fitting) the initial 
conditions for the hydrodynamical evolution of the sys- 
tem. This model is often referred to as the "final-state 
saturation" model since the gluon saturation here refers 
to the saturation in gluon production rather than to that 
in the initial state (color-glass condensate) wave func- 
tions of the colliding nuclei. The purpose of this paper is, 
in addition to the implementation of the next-to-leading- 
order (NLO) pQCD updates, to improve the conceptual 
and technical NLO framework of the EKRT modeling, to 
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make an effort to quantify the underlying uncertainties, 
and to study the predictive power of the approach. We 
show that after the improvements this fairly robust mod- 
eling leads to a very reasonable phenomenology for the 
particle multiplicities and pt spectra simultaneously at 
RHIC and the LHC. 



A. Original pQCD + saturation + hydrodynamics 
framework 

Let us first look back at the key elements as well as 
the key predictions in the original EKRT model setups 
[H El- The gluon saturation in central (b = 0) 
A+A collisions was implemented through an uncertainty- 
relation based geometrical saturation criterion, 

Naa(po, v 7 ^, AF = 1, b = 0) x = K sat nR 2 A , (1) 

Po 

where Naa is the number of produced minijets (to a good 
first approximation all gluons) [f| which fall into the mid- 
rapidity acceptance window AY = 1 with a few-GeV 
transverse momenta, pt > Po Aqcd- The factors tt/pq 
and ttR\ (with Ra the nuclear radius) account for the 
transverse area occupied by the production of each mini- 
jet and that available in a central A+A collision, and the 
proportionality constant if sa t was set to unity. In terms 
of collinear factorization, i.e. DGLAP-evolved nuclear 
parton distribution functions (nPDFs) fi/Ai an( i P er ~ 
turbativcly computable leading-order subcross-sections 
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d<7 lJ ' ->fei the minijet number can be computed as 2 

Naa(po, Ay,0) = T AA (0)a(N} AY , Po , (2) 

where 1^(0) = J d 2 s T A (s)T A (s) is the standard nu- 
clear overlap function, expressed in terms of the nuclear 
thickness functions T A , and the hard QCD cross section 
is (for the parton bookkeeping details, see Ref. 



v(N)ay, Po = Kqgb J dp' T dy 1 dy 2 S 2 (N;p 1 ,p 2 ) x 



x ^2x 1 f l/A (x 1 ,Q 2 )x 2 f J / A (x2,Q 2 ) 

ijkl 



da l ^ kl 
di 



(3) 



The " measurement function" S 2 [12|, |l3| above contains 
the imposed pr cut-off po and rapidity acceptance AY, 

S 2 (N; Pl ,p 2 ) = e(p T > po)[0(vi e Ay) + % 2 g AY)}, 

(4) 

where the labels 1 and 2 refer to the final state par- 
tons and pt to their (identical) transverse momentum. 
The functions denoted by 8 are the usual step functions. 
The factor Kqcd in Eq. (j3]) above accounts for the NLO 
corrections. In the original EKRT paper Q a constant 
Kqcd = 2 was assumed for simplicity (and since no cal- 
culation for these existed) but in the later setups [T3, El 
a y/sN jv-dependent .Kqcd was introduced on the basis of 
NLO computations of minijet Et production (l2l. Il3j. 

Once the saturation momentum p Q = p sat fulfilling the 
criterion in Eq. (TTJ is found, the amount of transverse 
energy carried by the minijets (whose px > Psat) into the 
rapidity acceptance Ay can be computed as 



E AA ( Po , 7^7, AY, 0) = T aa {0)<j(E t )ay, 



(5) 



where ct(Et) AY,p is obtained from Eq. ([3]) by replacing 
the measurement function S 2 (N;pi,p 2 ) with a new one 
for Et, 

S 2 (E T ; Pl ,p 2 ) = e( PT > P o)p T [0(yi € AY)+0{y 2 e AY)]. 

(6) 

The initial conditions for hydrodynamics are then ob- 
tained from the computed £t(po = p S nt) and satura- 
tion scale p sa t by assuming that the system thermal- 
izes essentially at formation, tq — l/p S at- The aver- 
age initial energy density becomes then (see [3, [HI [ll[ ) 
(e) = E AA / (-kR\tqAY). As discussed in [3], these initial 
conditions fix the rapidity density of entropy dS/drj in the 
system, which in turn, with ideal hydrodynamics, is di- 
rectly proportional to the final-state multiplicity. In the 
original EKRT framework [3] , where hydrodynamics with 
only 1-dimcnsional Bjorken scaling flow was considered, 
the main prediction was the obtained scaling law of the 



particle multiplicity, dN/dy oc A°- 92 s% 1 ^'" ' 20 , where the 
binary collision scaling A 4 / 3 of independent minijet pro- 
duction was thus tamed to be close to a wounded-nucleon 
scaling, and where a definite power-law behaviour in 
y/sNN was predicted. Given the robusteness of the model 
and also that these predictions were prepared before any 
RHIC data were available, the obtained ^/snn scaling 
turned out to be surprisingly close to the one obtained 
with the first LHC heavy-ion ALICE data 0]. Another 
interesting observation in the original EKRT model [3] 
was that the ratio E AA (p sat )/N AA (p sat ) is very close to 
the ratio e(T)/n(T) of an ideal massless boson gas, which 
suggests that further gluon multiplication is not neces- 
sary for (kinetic) thcrmalization. This observation, to- 
gether with the fact that the system becomes overdense 
with gluons at saturation, lends support to the early ini- 
tialization time 1 /psat of the hydrodynamical evolution. 
In the later, more detailed EKRT setups [ill EH 
where the pQCD calculation induced the ^/snn- 
dependent -Kqcd factors [T3, HH mentioned above, the 
saturated minijet initial conditions served as input for 
ideal 1+1 D hydrodynamics (central collisions, azimuthal 
symmetry, longitudinal scaling flow but dynamically 
evolving transverse flow) once a binary-collision (BC) 
transverse profile for the energy density, 



e(s) = T A (s)T A (s) 



a{E 



T)AY,p m 



r AY 



(7) 



was assumed. Also a more realistic equation of state, 
resonance decays and centrality selection were consid- 
ered in the hydrodynamical description. With such a 
setup, the charged-particlc multiplicities measured in 
central Au+Au collisions at RHIC at v /s]vjv = 56, 130 
and 200 GeV were all predicted quite nicely, see the 
comparison with the data in Refs. (To" 
ably, given that the prediction in Ref. 
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Remark- 
was made 

a decade before the LHC data and before the RHIC 
ytspjN — 200 GeV data, also the corresponding multi- 
plicity in central y^ijvw = 2.76 TeV Pb+Pb collisions at 
the LHC was predicted reasonably well. 3 

As discussed in Refs. [HI, [3], a good agreement with 
the measured pion and kaon transverse momentum spec- 
tra in central Au+Au collisions at RHIC (-^/sjvjv = 
130 and 200 GeV) was also obtained. The first EKRT- 
based results obtained in Ref. [l6| for charged- hadron pt 
spectra in Pb+Pb collisions at the LHC suggest that the 
framework can be expected to work well also there. What 
should be emphasized with these RHIC and LHC results 
is the robustness of the framework: due to the pressure, 
the PdV work, the rapidity density of energy degrades by 
almost a factor three during the hydrodynamical evolu- 
tion, so that it is quite nontrivial that a good agreement 



2 We adopt the notation of Refs. [TH . Il3tl here. 



3 The result obtained by interpolation from the predictions shown 
in Ref. [TTJ] agreed with the ALICE multiplicity (upper limit of 
the experimental error bar) within 7%. 
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with the measured px spectra follows. 

Also the centrality dependence of multiplicities in the 
EKRT approach has been studied, using the optical 
Glauber model for the collision geometry and localiz- 
ing the saturation criterion of Eq. ((T|) in the transverse- 
coordinate plane [I?)- As pointed out in Ref. [lH, a 
good agreement with the RHIC data for multiplicity-per- 
partipant-pair vs. the number of participants is found 
when the same (optical) Glauber model is applied also 
in the data analysis - for comparison, see Fig. 23(a) in 
Ref. [H] and Fig. 4 in Ref. [ljj ■ Elliptic flow predictions 
in the EKRT framework can be found in Ref. [2(| for 
RHIC and in Refs. [HS2 for the LHC. 



B. Open questions with the original EKRT setup 

Inspite of the working phenomenology achieved, there 
are a number of open questions and shortcomings with 
the original EKRT model setup which we will consider in 
the present study: 

(i) The formulation of the model with saturation of the 
number of produced gluons, Naa in Eq. (fTJ), is problem- 
atic, since the number of gluons cannot be defined in a 
manner which would be collinear (CL) and infra-red (IR) 
singularity safe also in NLO pQCD: without introducing 
an extra resolution scale, how should one count soft glu- 
ons and two collinear gluons? Consequently, the factors 
Kqcd in Eq. ([3]) cannot be computed directly for Naa 
but they have been extracted from an NLO calculation 
of transverse energy E AA which is a CL/IR-safe quan- 
tity and can thus be rigorously defined and computed 
as in Refs. [lU, EH- These problems would obviously be 
avoided if the saturation would be required for E AA in- 
stead of Naa- 

(ii) There are uncertainties related to the geometri- 
cal final-state saturation criterion: It has not been clear 
whether an explicit strong coupling constant a s should 
appear in Eq. (TTJ) if it describes a fusion of the produced 
gluons. Also, it has been unclear what is the role of the 
possible further proportionality constant K sat in Eq. ((T|) . 
Also, if the softer gluons at px < Psat are produced but 
they fuse in the overdenscly populated phase space, one 
may ask whether the total transverse energy should still 
increase with the produced soft gluons? 

(Hi) As explained above, the hydrodynamical evolution 
in this framework is initiated with the produced minijet 
transverse energy and not with their number, computed 
at pa = p sat once p sa t has been determined on the basis of 
the minijet number. The whole procedure of computing 
the produced initial energy densities for hydrodynamics 
would obviously become much more straightforward if 
the saturation criterion could be formulated for the pro- 
duced E AA instead of Naa- 

(iv) The saturation criterion in Eq. (p} is extensive in 
AF on the left-hand side but not on the right-hand side 
(see also the discussion in Ref. [22| ). Since the previous 
studies have been made for AY = 1, this has not been a 



problem in practice but obviously also this question calls 
for an improved formulation of the saturation criterion, 
where Ay would explicitly appear also on the right-hand 
side of Eq. (HJ. 

(v) The earlier EKRT studies (as well as the NLO E T 
studies [H) have so far applied only the EKS98 [H[ LO 
nPDFs in the computation of the minijet cross sections. 
An NLO update regarding the nPDFs should be done, 
and also a study of the propagation of the nPDF un- 
certainties into the hydrodynamic initial conditions and 
final state multiplicities should be performed. 

In this paper, we address the above issues as follows: 
First, we reformulate the gluon saturation in terms of 
transverse energy per rapidity unit. This solves the first 
item above by removing the need for the phenomeno- 
logical factors -?Cqcd for Naa, since we can now directly 
compute the saturating quantity in NLO pQCD. Also the 
third item is then solved. To address the second item, we 
show how the geometric saturation criterion without an 
explicit a s (or its powers) arises with a more detailed de- 
scription of the EKRT saturation mechanism. Also the 
fourth item gets conveniently solved with this reformula- 
tion. To address the fifth issue, we bring the NLO com- 
putation for minijet Et pro duction up-to-date by using 
the EPS09 NLO nPDFs [H. We also study the nPDF- 
originating uncertainties in the final-state multiplicities 
using the EPS09 error sets. 

Our aim here is to address also the uncertainties and 
phenomenological parameters of the approach more con- 
cretely than before. In the rigorously computable pQCD 
part, we discuss the freedom in defining the measure- 
ment function(s) for the NLO Et calculation. In the 
more phenomenological saturation part, we quantify the 
uncertainty due to the unknown proportionality constant 
in the saturation criterion. We will also show how these 
uncertainties (or rather, freedoms) are correlated when 
studying the final-state multiplicities at the LHC and 
RHIC. At the same time, we chart the predictive power 
of the improved EKRT framework. 

Once getting the pQCD+saturation part and multi- 
plicities under improved control, we also consider the 
bulk hadron px spectra in central Au+Au collisions at 
RHIC and Pb+Pb at the LHC. To get an updated EKRT 
baseline for further improvements, we apply ideal hydro- 
dynamics with a statc-of-the art equation of state. It 
should be emphasized that our goal here is not to tune 
the model to fit the px spectra at RHIC and LHC as 
perfectly as possible but, rather, to study whether a rea- 
sonable agreement especially with the LHC px spectra 
[25| can be found, by keeping exactly the same parameter 
setup for RHIC and LHC. We also discuss how the im- 
proved understanding of the EKRT framework presented 
in this paper leaves room also for further improvements 
such as adding dissipation (viscosity) into the hydrody- 
namical description. 

The rest of the paper is organized as follows: In Sec. [H] 
we present the improvements and updates of the NLO 
pQCD + saturation framework and specify our hydrody- 
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namical setup. Section IlIII contains the obtained results 
for the systematics of multiplicities and pr spectra at 
RHIC and the LHC. In Sec. IIVI we summarize and dis- 
cuss the further improvements left for future work. 



II. IMPROVED EKRT FRAMEWORK 
A. Saturation in Et 

We take the following new angle in interpreting the sat- 
uration in the EKRT framework: instead of a saturation 
of the number of produced final state gluons we suggest 
the saturation to take place in transverse energy produc- 
tion when the 3 — > 2 and higher-order (n > 4) —> 2 pro- 
cesses start to dominate over the conventional 2 — > 2 pro- 
cesses. The saturation should then not be considered just 
as final (or initial) state saturation but all (n > 3) — > 2 
processes which decrease the gluon number relative to 
the independent 2 — > 2 processes, are in effect. Also, as 
all the higher order processes then effectively reduce the 
number of final-state gluons, and since the initial state 
gluons do not carry transverse momentum, the produc- 
tion of Et from the collisions at pt < p S at is conjectured 
to be negligible. 

We thus require that at saturation (which is assumed 
transversally non-local here), the rapidity densities of the 
transverse energy fulfill the condition 



-(2->2)~ -7^(3 



dy 



dy 



(8) 



Considering only the leading-order order cases in a s , we 
assign a factor Ta • g for each of the incoming gluons 
(g for gluon PDFs, and T A ~ A/(-kR 2 a )), the factor nR\ 
for the transverse integration d 2 s in Taa, the appropriate 
powers of a s , the pr cut-off scale po for the Et, the scale 
Pq 2 for <t(2 — s- 2), and the scale p^ 2 to compensate for 
the fm~ 2 dimension of the extra Ta in the 3 —J- 2 case, 
we arrive at a scaling law 



TrR 2 A (T A g) 2 ^Po ~ irR 2 A (T A g) 3 ^Po- (9) 
Po 



l_a| 

2 2 - 

Po Po 



At saturation, po = Psatj this leads to a scaling 



Ta9 



Pi 



(10) 



for the gluon density probed at saturation. 4 Feeding this 
scaling back to Eq. ([5]) gives 



dE 7 
dy 



(2 -> 2) cx R\pl 



(11) 



4 Interestingly, a similar relation is traditionally obtained in the 
CGC framework, see e.g. Ref. [2^1 . 



at saturation. Thus, we arrive at the following 
geometrical-like saturation criterion for the average mini- 
jet transverse energy produced in a central A+ A collision: 

E$ A (p , 0) = K sat R 2 A p 3 AY, (12) 

where no explicit a s appears, the rapidity interval AY 
(from dy) appears also on the right-hand side, and where 
the proportionality constant A" sat is to be determined on 
the basis of the measured data at one chosen cms-energy 
This is the saturation conjecture we now test in the fol- 
lowing. 



B. NLO computation of minijet Et 

The exact NLO formulation of minijet Et produc- 
tion, introduced originally in Refs. [Ill, [l3|, is based on 
collinear factorization and the subtraction method [27j . 
As shown in Eq. (|5|), the evaluation of the produced Et 
in central A+A collisions is carried out by computing 
the first Et moment <j{Et) AY,p of the perturbative Et 
distribution of minijets in a single N+N collision and in- 
cluding the nuclear collision geometry through the stan- 
dard nuclear overlap function Taa(^>)- The NLO pQCD 
framework is described in detail in Ref. [l3[ but we briefly 
recapitulate the formulation of <j(Et) here, to discuss a 
generalization of the measurement functions for Et- 

The semi-inclusive Et distribution of minijets in a ra- 
pidity interval Ay in N+N collisions [!, [TU, [l3[ can be 
computed to NLO pQCD as 



da 



dE T 



da 



dE T 



+ 



(la 



2->2 



dE T 



AY, Po 



3 .mm i>3; 



d[PSh 



(13) 



where the integrations take place in 4 — 2e spacetime di- 
mensions, and the 2 — > 2 and 2 — > 3 differential partonic 
cross sections are denoted as 



da 



2^2 



da 



2-S-2 



d[PS} 2 dp T 2dyidy 2 d(f)2 
da 2 ^ 3 da 2 ^ 3 



(14) 



d[PS] 3 dyidy2dy 3 dpT2dpT3d(/>2d(j)3 



For the two-parton final state the appropriate kinemati- 
cal variables are the rapidities yi, j/2, transverse momen- 
tum pT2 and azimuth angle <f>2 ■ From transverse momen- 
tum conservation, we have pti = PT2 and <p± = 4>2 + 7r 
as we do not include any intrinsic transverse momen- 
tum. Similarly, for the three-parton final state, the suit- 
able kinematical variables are yi,y 2 ,V3,PT2,PT3, </>2 and 
03, since the transverse momentum conservation fixes 
Pti = — (PT2 + PT3)- The measurement functions 
S 2 = S 2 {pi 1 p 2 ) and S 3 = S 3 (pi,p 2 ,P3), which depend 
on the four-momenta of the final-state partons, define 
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the physical quantity to be computed. In this case this is 
the Et distribution of minijets which fall into a given ra- 
pidity acceptance window Ay and which originate from 
hard (pcrturbative) collisions where at least an amount 
2p of transverse momentum is produced. Thus, the mea- 
surement functions also define what we mean by a hard 
process here. 

To regularize the IR/CL divergencies present in the 
partonic NLO cross sections, we must consider the 
ultraviolet (UV)-renormalized squared 2 — > 2 and 2—^3 
scattering matrix elements of order af, in 4 — 2e dimen- 
sions and MS scheme. The divergent terms show an e _1 
and e~ 2 behavior. The full analytical calculation for 
these matrix elements was done in Rcf. 28] 5 . The ex- 
plicit cancellation of the IR/CL divergencies takes place 
only if the three-parton measurement function, S3, re- 
duces to the two-parton one, S2, in the IR and CL limits 
[27| , i.e. when one of the final state partons becomes soft, 
or collinear with any other parton in the process. 

In a hard scattering of partons to NLO, we may have 
one, two, three or zero minijets in the rapidity acceptance 
region Ay, which here is the mid-rapidity unit, 



Ay : \y\ < 0.5, < cj> < 2ir. 



(15) 



All the partons are assumed massless, thus the transverse 
energy within Ay is the sum of the absolute values pTi of 
the transverse momenta of those partons whose rapidities 
are in Ay: 



where the step function e(yi) is defined as 



1 ify.eAy 
otherwise. 



(16) 



(17) 



In the LO and NLO 2 — > 2 cases, where the transverse 
momenta arc equal in magnitude, pti = PT2, the hard 
scatterings can be defined to be those with large enough 
transverse momentum, regardless of where the partons 
go in rapidity: pt > Po S> Aqcd, or equivalently, 



Pti + Pti > 2p , 



(18) 



where po is a fixed external parameter which docs not 
depend on Ay. This readily generalizes to the NLO 2 — > 
3 processes as 



Pti + PT2 + PT3 > 2po- 



(19) 



As discussed in Ref. [l3| , a possible further element in 
defining the measurement functions here is that in the 
2 — > 3 case we may still restrict the amount of minimum 
Et at Ay in an IR/CL safe way: In the 2 — > 2 case, the 



5 Details of some of these rather complicated calculations will be 
elucidated in 1291. 



non-zero Et in Ay is always larger than po, while in the 
2—^3 case we can have (and have plenty of, see [l3j]) hard 
processes where two partons fall outside Ay and one soft 
parton inside. At the IR limit in this special case, we 
obviously have no Et in Ay and the usual 2 — 2 limit 
is correctly recovered. The other equally well IR/CL safe 
extreme case is that we could require the Et in Ay be 
at least po as always is in the 2 — > 2 case. The new 
feature introduced in the present study is that in fact any 
minimum amount of Et between and po constitutes an 
equally good, IR/CL safe restriction for the Et in Ay 
which relaxes back to the 2 — > 2 case at the IR and CL 
limits. 

By combining the definitions of Et in Ay in Eqs. (|16[) . 
(TT8"|) and (jT9")l together with the definition of the hard 
scatterings and the possible restriction of Et discussed 
above, the IR/CL safe measurement function S2 can now 
be written down for the 2—^2 scatterings as 



(20) 



S2(Pl,P2) = 0(PT1 + PT2 > 2p ) 

5(E T - {e(yi)pTi + e(2/2)PT2}) 
and 5*3 for the 2 — >■ 3 scatterings as 

Sd,{pi,P2,P3.) = 0(pti +PT2 +PT3 > 2p )6(E T > /3p ) 
5(E T - {t{yi)pTi + e(y 2 )pT2 + e(?/3)PT3}), 

(21) 

where the constant j3 € [0, 1]. The measurement function 
53 above thus generalizes the formulation of Ref. [l3| . 
where only the special cases f3 ~ and 1 were studied. 
In the IR and CL limits, we can see that S3 — > S2 as is 
required for the cancellation of the corresponding diver- 
gencies. Importantly, we notice already here that with 
the parameter (3 we can control how much Et is allowed 
to form within a certain formation time oc 1/poi which in 
turn affects the hydrodynamical initial energy densities. 
Thus, (3 is a parameter of this model, whose value is to 
be determined from the RHIC and/or LHC data as will 
be explained in Sec. III. 

From Eq. (|13p we now obtain the first moment of the 
semi-inclusive Et distribution as 



a(E T ) 



Po.AY 



(LEt Et 

dE T 



PoAY 



(22) 



TIpoAY 



■a(E> 



TIpoAY 



where, in 



= 777 / d l PS h 



da 2 



d[PS\ 



-S 2 {pi,P2) (23) 



and 



->3 
,AY 



d[PS}. 



do~ 



2^3 



d[PS] 



-S 3 (pi,P2,P 3 ), (24) 



we have integrated the delta functions away, so that our 
final measurement functions can be written as 



£2(^1, P2) = 



e (2/i) + e(2/2) 



PT2d{pT2 > Po) (25) 
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and 

S3(pi,P2,P3) = (e(yi)PTi + AV2)VT1 + e(y 3 )p T 3) 

d{PTl + PT2 + PT3 > 2p ) 

Q{t{Vl)PTl + t{V2)PT2 + ^{yz)PT3 > Ppo)- 

(26) 

Naturally, also the measurement functions 5*2 and S3 ful- 
fil the IR/CL-safety criteria and thus the Et moment 
c(Et)po,ay is a well-defined IR/CL safe quantity to com- 
pute. 

As a straightforward improvement of the EKRT frame- 
work, we now employ the NLO EPS09 nPDFs [H| and 
CTEQ6M [H PDFs. In fixing the renormalization scale 
Hn and factorization scale /if, we follow the common 
practise and choose them to be equal, pr = /if = /i- We 
set fj. to be proportional to the hardness of the collision, 
i.e. to the total transverse momentum produced in the 
hard process, regardless of the partons being in Ay or 
not: 



/' 
/' 



N^PTl 
N^pri 



■p T2 )/2 = N^pt 

PT2 +PT3)/2, 



(27) 



where we choose = 1. Note that this choice of /j is 
IR/CL safe, as is required for the exact cancellation of 
the divergences. 

Thus, for the very first time, we are now able to com- 
pute the minijet Et production in A+A collisions as rig- 
orously as currently possible to NLO pQCD, and supple- 
ment this calculation with the saturation of Et as de- 
scribed above. We will also study how the nPDF uncer- 
tainties, described by the error sets in EPS09, propagate 
into the computed Et, and thus to the hydrodynamical 
initial conditions. 



C. Hydrodynamical setup 

In this work we use ideal hydrodynamics to evolve the 
QCD-matter initial state given by the EKRT model to 
final state particles. We solve the equations 



0. 



(28) 



where = (e + P)u^u" - Pg^ v is the energy- 

momentum tensor, e is energy density, P is pressure and 
is the fluid four-velocity. In this work we are only 
interested in the multiplicity and pj-spectra in central 
A+A collisions and thus we can make some approxima- 
tions in order to speed up the calculations. Firstly, we 
assume longitudinal boost-invariance which is a good ap- 
proximation when discussing results for mid-rapidities. 
Secondly, since we limit our studies to central collisions, 
we can assume azimuthal symmetry and thus the original 
3+1-dimensional numerical problem is reduced to 1+1 
dimensions. In addition, we have set the baryochemical 
potential to zero. We employ the SHASTA algorithm 
[3l| for solving the hydrodynamical equations. We use 



s95p-PCE-vl equation of state (EoS) [32j (PCE for par- 
tial chemical equilibrium) to close the set of equations. 
This EoS has a chemical freeze-out at T c h om = 150 MeV. 

During the hydrodynamical evolution we construct an 
isothermal freeze-out hypcrsurface E M using a criterion 
Tf = 120 MeV, which was chosen so that we get a reason- 
able agreement with the measured pt spectra at RHIC. 
We assume that the kinetic freeze-out happens instan- 
taneously on this surface and the fluid is converted to 
particles using the Cooper-Frye formula \-VA . After we 
have obtained the thermal particle spectra, we sample 
particle ensembles from them in the same way as was 
done in Ref. 34[. The rapidity of the particle is taken 
from a flat rapidity distribution using an interval \y\ < 3. 
Strong and electromagnetic 2- and 3-particle decays are 
then done one particle at a time. 

Since the transversally-averagcd EKRT model consid- 
ered here does not fix the transverse profile for the pro- 
duced initial energy density, we have employed different 
profiles to get a handle on the uncertainties related to 
the profile. Our choices are binary collision and wounded 
nucleon (WN) profiles from the optical Glauber model. 
With the BC profile, the initial state is obtained as in 
Eq. ([7]) and with the WN profile the initial energy den- 
sity is 



(29) 



e(s) =K W n[T a (s)(1 - exp[-a NN T B (s)}) 
+ T B (s)(l - cxp[~a NN T A (s)]fj , 



where the overall normalization constant -Kwn is fixed 
so that we have the same amount of energy, ETijps&t), 
as with the BC profile. Since the entropy is obtained by 
converting the energy density to entropy density using 
the EoS, the total amount of entropy, and thus also the 
final multiplicity, is slightly different with the WN profile 
than with the BC profile. 

The hydrodynamic initial time is always taken to be 
To = 1/psat- This relation contains a further 0(1) pro- 
portionality constant which we set to unity throughout 
this work. For the phenomenology, the value of this con- 
stant is obviously correlated with the parameters K sat 
and (3 but a quantitative investigation of this is left for 
future studies. 



III. RESULTS 

Next we discuss the results from the present NLO- 
improved pQCD + saturation + hydrodynamics frame- 
work. We start with the computed NLO minijet Et, 
charting in particular the effects of the free parameters 
/3 and K sat , as well as the NLO nPDF uncertainties. We 
demonstrate how the i/p-saturation works at RHIC and 
LHC, and what is the effect of the above uncertainties 
on saturation, i.e. on the hydrodynamic initial condi- 
tions. After understanding these, we fix a specific (/3, 
Ksat) combination, and compare our results with the 
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LHC and RHIC data on charged-particlc multiplicities 
as well as on identificd-hadron pt spectra. The propa- 
gation of the nPDF errors into our multiplicity results 
is charted in detail. As will be seen, the outcome for 
the predictive power of the improved EKRT approach is 
quite encouraging. 

Figure [TJ shows the average NLO minijet transverse 
energy E AA produced in 0-5% central A+A collisions in 
the rapidity acceptance Ay = 1 as a function of the 
cut-off po- We consider Au+Au collisions at the RHIC 
energy ^snn = 200 GeV and Pb+Pb collisions at the 
LHC energy y/sjyN = 2.76 TeV. The centrality selection 
is simulated by considering a central A e s+A e s collision 
of an effective nucleus, A c g = 181 at RHIC and A e ff = 
193 at the LHC, as explained in detail in Ref. The 
NLO E AA curves are computed from Eqs. ([5]) and (|22|) 
using the measurement functions (|2"5|) and ([2"rJ]) with the 
parameter (3 = 0, 0.5, 0.75 and 0.9. As seen in the figure, 
an increasing f3, which cuts the accepted phase space via 
the restriction Et > /3po, can reduce the amount of the 
produced Et as much as by a factor ~ 2. 

The rising curves in Fig. [TJ are the right-hand side of 
the saturation equation (fT2"j) . K sat R\ c!! pQAY , with the 
proportionality constant K sat = 0.5, 1, and 1.5. 6 Sat- 
uration thus takes place when these curves cross those 
for Et- As the figure demonstrates, the values of j3 and 
-Ksat are correlated: the same Et can be obtained with 
many different K sa ±) pairs. External input - experi- 
mental data - is needed to resolve the best allowed val- 
ues for these parameters. To get a hold on this, the red 
bands in Fig. Q]are to indicate which E^: A (po) would be 
needed according to our hydrodynamic prescription, as- 
suming To = 1/po and the BC energy-density profile, to 
reproduce the experimentally measured multiplicities at 
RHIC and LHC. 7 The solid lines inside the red bands are 
computed from the statistically weighted averages of the 
measured charged-particle multiplicities at the LHC (AL- 
ICE [3 and CMS and RHIC (PHENIX [H, STAR 
H3 and BRAHMS [38j|). The red error bands correspond 
to the largest (smallest) values of the experimental up- 
per (lower) limits of multiplicities, mapped now into the 
QCD-matter initial state through ideal hydrodynamics. 

We can also see from the figure that if we choose f3 = 
and tune K sat to fit the RHIC data, we overshoot the " ex- 
perimental" error band at the LHC - this was a problem 
with the old EKRT setup, where the yG scaling of the 
multiplicity became too strong. Interestingly, however, 
we observe that one effect of our hardness-parameter j3 
is the taming of this yfs scaling: if we choose, as perhaps 
the most natural choice would be, K sat = 1 (the solid ris- 
ing curve) and (3 = 0.75 (the solid black Et curves), we 



6 The A e ff is different for RHIC and LHC, hence the splitting of 
these curves in the figure. 

7 Thus, the input for the red bands contains only hydrodynamics 
and experimental data, no pQCD or saturation. 



match the average LHC multiplicity perfectly and agree 
very nicely also with the RHIC average multiplicity. As 
seen in the figure, the acceptable values of (3 and K &B ± 
are clearly correlated: with a larger (3 a larger value of 
K sat is needed. 
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FIG. 1: (color online) The NLO minijet transverse energy 
E AA produced in 0-5% central A+^4 collisions in the mid- 
rapidity region AY = 1, as a function of the pr cut-off po- 
Upper set of the Et curves: LHC Pb+Pb at ^/sjvjv = 2.76 
TeV with p = 0, 0.5, 0.75 and 0.9. Lower set: RHIC Au+Au 
at = 200 GeV with the same values of j3. The rising 

curves arc K Ba ,tR 2 A e[f PoAY, the r.h.s. of Eq. (TJJ}, with if 8at = 
0.5, 1 and 1.5, and°yl ff = 193 (181) for the LHC (RHIC). The 
red bands labelled as "Data+Hydro" show how the measured 
multiplicities translate into the initial states according to ideal 
hydrodynamics. 

Figure [2] shows the computed NLO E^ A in Ay = 1 for 
fixed (3 = 0.75 at RHIC (lower bands) and LHC (upper 
bands). The straight lines in this log- log plot are again 
the r.h.s. of Eq. (|12[) with a fixed K sat = 1, and the A e gs 
simulating the 0-5% ccntralities at RHIC and LHC. The 
black lines inside the error bands are the E A computed 
by using the EPS09 best fit S°. The yellow error bands 
show the uncertainty which originates from the nPDF 
uncertainties, computed using the 30 error sets in EPS09 
according to the prescription given in (24|: 



(AE AA )± = (max(£^)±) 



where 

max(^)+ 



lA (S+) - E AA (S°) 
E AA (S k 



max{E T 

^ AA {S t 



and 



max(E AA 



\- =max{E AA (S°) 
E AA (S°) 



E AA (S°),0} 

■^ A (s+), 

E AA (S^),0}, 



(30) 



(31) 



(32) 
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FIG. 2: (color online) The computed NLO minijet E? A 
in AK = 1 with fixed /3 = 0.75, for RHIC Au+Au colli- 
sions at y'SNN — 200 GeV and LHC Pb+Pb collisions at 
^snn = 2.76 TeV, including the nPDF uncertainties from 
the 30 EPS09 error sets. The ArffS simulate the 5% most 
central collisions. 



tainties are smaller than in the smaller-x region probed 
at the LHC. 

Finally, and perhaps most importantly, Fig. [5] shows 
that the rigorous NLO Et computation itself seems to 
be well under control even at the rather low perturbative 
scales at saturation, and, that the nPDF-originating un- 
certainties remain rather modest, ca. 15 % at saturation 
both at the LHC and RHIC (see Table 1). 
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where E AA (S^) denotes the value of the E AA computed 
with the set Eft for k = 1, ... , 15. 

In Fig. [21 we observe that at the LHC the error bands 
shrink consistently towards higher po . This is due to the 
DGLAP scale evolution of nuclear gluon PDFs, which 
rapidly decreases the shadowing of small-x gluons as well 
as their uncertainties. Also note here that the minimum 
scale of both the CTEQ6M PDFs and EPS09 nuclear 
effects is Qo = 1.3 GeV, and that for the nPDFs below 
Qo we have used those at Qo for simplicity. Thus the 
results at po < 1-3 GeV in Figs. [1] and [2] are not fully 
consistent with the DGLAP evolution. However, with 
the choice /3 = 0.75 and K sat = 1, we conveniently have 
Psat > Qo at saturation also at RHIC (sec Table 1 ahead). 

For the RHIC results in Fig. [5J we notice a more non- 
monotonic behaviour in the widths of the corresponding 
error bands, as well as a few-percent numerical uncer- 
tainty at 1.3 < po < 2 GeV. These can be traced back to 
the non-trivial behaviour of the large NLO 2 — > 2 contri- 
bution and its interplay between the 2—^3 contributions 
which require 6 dimensional MC integrations [l3l |. Sim- 
ilarly, the small numerical fluctuations seen at the LHC 
curves originate from these multidimensional MC inte- 
grations. 

Furthermore, we can see a general trend that for a fixed 
Po > 2 GeV the error bands are smaller at RHIC than at 
the LHC. This again is a reflection of the different mag- 
nitudes of the uncertainties in the nuclear gluon PDFs at 
different values of x: at RHIC, due to the smaller ^/snn, 
one probes larger values of x where the shadowing uncer- 



FIG. 3: (color online) Charged-particle multiplicity dN c h/drj 
at \r)\ < 0.5 in 0-5% most central y/sJ^N = 2.76 TeV Pb+Pb 
collisions at the LHC and ^/snn = 200 GeV Au+Au collisions 
at RHIC, obtained from NLO pQCD + saturation + hydro- 
dynamics using the BC and WN initial energy-density profiles 
(filled red and blue squares, correspondingly), and a compar- 
ison with the ALICE and CMS data (LHC) Refs. [1 |al and 
PHENIX, STAR and BRAHMS data (RHIC) Refs. JM^H]. 
Notice that most of the points have been shifted sideways to 
improve their visibility. 

In Fig. [3] we show the computed charged-particle mul- 
tiplicity dNch/dr) at \r]\ < 0.5 in 0-5% most central 
y/sW^ = 2.76 TeV Pb+Pb collisions at the LHC and 
y/spjN = 200 GeV Au+Au collisions at RHIC, obtained 
in the current framework of NLO pQCD + saturation + 
hydrodynamics with /3 = 0.75 and K sa ,% = 1, using the 
BC and WN initial energy-density profiles. Also the com- 
parison with the ALICE [land CMS H datajLHC) 
and PHENIX STAR [37| and BRAHMS [13 data 
(RHIC) is shown. 

The error bars to our results in Fig. [3] have been ob- 
tained as follows: First, we calculate the NLO E AA in 
AF = 1 for RHIC and the LHC by fixing (3 = 0.75 
and using the best fit and the 30 error sets of EPS09. 
Second, we solve the saturation equation (T2"|) and deter- 
mine p sa t with a fixed K sat = 1 for all the 31 different 
NLO E AA results, and thus at saturation we find 31 dif- 
ferent (-Ey^PsatjjPsat) pairs. Third, we construct the 
hydrodynamical initial conditions, e(s, To, & (Et) p aa , t ) at 
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r = 1/psat (see Eqs. ([7]) and (j2"9"j) ) and calculate dN^/dr] 
for every EPS09 set separately. Finally, using these num- 
bers, we compute the theoretical error bars using the 
EPS09 prescription of Eq. ([30jl . We can see that the 
error bars to our results are slightly larger than the ex- 
perimental error bars at the LHC, and of the same order 
as those at RHIC, indicating again that the calculation 
presented here is not suffering from a large gluon nPDF 
uncertainty. 

Thus, from Fig. [3] we conclude that our NLO pQCD 
+ saturation + (ideal)hydrodynamics framework repro- 
duces the measured RHIC and LHC multiplicities quite 
nicely and the ^snn scaling seems to work very well. 
The transverse-profile uncertainty in the computed mul- 
tiplicity is a few percent and the nPDF-related uncer- 
tainty about ±10%, for fixed (5 and K sat , both at RHIC 
and LHC. Most importantly, since it is possible to keep 
the values of /? and K sat fixed, we can conclude that our 
framework has some definite predictive power. 
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TABLE I: The collision parameters, the obtained NLO pQCD 
+ saturation key results for fixed /3 = 0.75 and K sa _ t = 1, the 
corresponding hydrodynamical input and final results for the 
charged-particle multiplicities. 



Tabic 1 collects the collision parameters as well as 
our results for p sat , -E^^Psat) at saturation, the corre- 
sponding hydrodynamic inputs To and ct(Et}, as well as 
the multiplicities. The errors shown for these quantities 
(upper and lower limits separately) have been computed 
based on the 30 error sets of EPS09. 

Figure 2] shows the pt spectra of tt + , K + , p and p in 
0-5 % most central Au+Au collisions at ^snn = 200 
GeV at RHIC, computed from NLO pQCD + satura- 
tion + hydrodynamics with (3 = 0.75 and -Ksat = 1, and 
the WN and BC initial energy- density profiles. These 
results thus correspond to the computed RHIC charged- 
particle multiplicities shown in Fig. |3l and we see that 
with the adopted PCE EoS and T f = 120 MeV, both 
the particle multiplicities and their px spectra measured 
by PHENIX STAR go[ and BRAHMS H3 become 
reproduced quite well. Note, however, that since we do 
not include any net-baryon number in the hydrodynami- 
cal calculation, wc cannot address the difference between 
the protons and antiprotons here. 
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FIG. 4: (color online) The pr spectra of tt + , K + , p and p in 
0-5 % most central Au+Au collisions at ^/snn = 200 GeV 
at RHIC, computed from NLO pQCD + saturation + hy- 
drodynamics with the WN profile (solid blue lines) and BC 
profile (dashed red lines), and measured by the PHENIX 39] 
(squares), STAR [U (circles) and BRAHMS (triangles) glf 
experiments. The calculated spectra and the PHENIX data 
are without the hyperon feed-down contributions. 



Figure [5] shows the corresponding identified particle 
Pt spectra for 0-5 % most central Pb+Pb collisions at 
y/SNN = 2.76 TcV, computed in our improved NLO 
EKRT framework with the same T f = 120 MeV as at 
RHIC, compared with the very recently published spec- 
tra from the ALICE experiment [42]]. Again, the com- 
puted spectra correspond to the computed LHC charged- 
particle multiplicities shown in Fig. [3] Very interestingly, 
we find that also at the LHC we can reproduce both the 
particle multiplicities and their px spectra quite well. Es- 
pecially, wc would like to emphasize that the computed 
results arc not a fit to these ALICE data but a predic- 
tion in the following way: First, we fixed the parameters 
j3 = 0.75 and A' sat = 1.0 by requiring that the measured 
charged-particle multiplicity is reproduced at the LHC 
(Fig. [T] the red band). Simultaneously, also the RHIC 
multiplicity comes out well. Then we fixed the kinetic 
decoupling temperature Tf = 120 MeV on the basis of 
Fig. QJ and kept it fixed when computing the LHC spec- 
tra in Fig. O 8 Regarding the initial energy-density pro- 
file uncertainty, we can conclude based on Figs. 0] and [5] 
that a flatter WN profile would seem to reproduce the 
identified particle px spectra better than the steeper BC 
profile. 



The charged-particle multiplicities are not sensitive to Tf. 
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FIG. 5: (color online) The pr spectra of tt + , K , p and p in 
0-5 % most central Pb+Pb collisions at ^/sjvjv = 2.76 TeV 
at the LHC, computed from NLO pQCD + saturation + hy- 
drodynamics with the WN profile (solid blue lines) and BC 
profile (dashed red lines), and measured by the ALICE [i^] 
(squares) experiment. 



IV. DISCUSSION AND OUTLOOK 

In this work, we introduced an improved formulation 
of the pQCD + saturation (EKRT) model Q, where 
the saturation of the gluon minijets occurs for the trans- 
verse energy, Et, rather than for the number of gluons. 
Since the minijet Et is a quantity safe from the infrared 
and collinear singularities, it can be computed rigorously 
to NLO pQCD as a function of the minijet transverse- 
momentum cut-off pq. The IR/CL safe minijet Et in 
a chosen acceptance region is defined by measurement 
functions analogous to jet physics. We showed that in 
fact a set of infinitely many equally possible IR/CL-safe 
measurement functions for the Et can be formulated in 
terms of a hardness parameter (3. In computing the mini- 
jet Et to NLO with these new measurement functions, 
we also updated the framework of Refs. [HI, HH with the 
EPS09 NLO nPDFs [Hj], and studied the propagation of 
the nPDF uncertainties to the minijet Et with the 30 
error sets of EPS09. Thus, in the present study the com- 
putation of the minijet Et production is performed for 
the very first time genuinely and consistently to NLO in 
A+A collisions. 

In the more phcnomcnological part of the study, we for- 
mulated a dynamical saturation criterion for the minijet 
Et , which solves the open issues of the old EKRT- model 
regarding the appearance of the rapidity acceptance in- 
terval AY" and powers of a s in the saturation criterion. 
The outcome is the EKRT-like geometrical saturation cri- 
terion for the collincarly factorized minijet Et produc- 
tion, containing one unknown proportionality constant 



.Ksat . Applying this saturation criterion to our NLO com- 
putation for the average minijet Et production in 5 % 
most central A+A collisions at RHIC and LHC, we de- 
termined the saturation scales p sa t and the amounts of 
Et produced into a mid-rapidity unit at saturation. In 
particular, we demonstrated how these depend on the 
model parameters /? and K saii , and we also quantifed the 
EPS09-originating error bars to these quantities. 

Converting the saturated minijet Et into the QCD 
matter energy density at tq = l/p sa t, thus assuming 
thermalization at formation, and adopting either a BC 
or a WN transverse profile, we obtained initial condi- 
tions for the hydrodynamic evolution. In this baseline 
study, where our main goal is to chart the general fea- 
tures and predictive power of the NLO-improvcd pQCD+ 
saturation + hydrodynamics modeling rather than a de- 
tailed fitting of the data, we restricted ourselves to boost- 
invariant ideal hydrodynamics with a state-of-the art 
PCE EoS [32 |. What we believe is particularly useful 
for understanding the initial-state phenomenology here, 
is that we concretely showed how the measured charged- 
particlc multiplicities translate into the QCD-matter ini- 
tial conditions through hydrodynamics (the red bands in 
Fig- Hi- Using this mapping, we could directly see which 
parameter pairs ((3,K slL t) - if any - would reproduce the 
measured LHC and RHIC charged-particle multiplici- 
ties simultaneously. Remarkably, the outcome is that 
(at least in this ideal-hydrodynamics framework) such 
a simultaneous reproduction of multiplicities is very well 
possible, and that we can determine the acceptable range 
of values for these correlated parameters. Using then one 
possible parameter combination, /3 = 0.75 and K sat = 1, 
we computed the identified particle pt spectra for the 5 
% most central A+A collisions at RHIC and LHC, fixing 
the kinetic decoupling temperature from the measured 
RHIC spectra. Again, the outcome is that the particle 
multiplicities and pt spectra are very nicely reproduced 
simultaneously both at RHIC and the LHC — empha- 
sizing the fact that the identified hadron pt spectra for 
the LHC Pb+Pb collisions which wc obtained here, is a 
prediction, and not a fit to the ALICE pt spectra. 

To summarize, we have shown that the NLO-improved 
pQCD + saturation + ideal-hydrodynamics is a viable 
model for describing particle production in central heavy- 
ion collisions at the LHC and RHIC. We have quantified 
the key parameters of the approach and studied the prop- 
agation of different uncertainties into the hadron multi- 
plicities and pt spectra. Most importantly, our results in- 
dicate that the framework has definite predictive power: 
The key-parameters of the NLO pQCD calculation (/?) 
and saturation (i^ sat ) can be fixed - not uniquely but in 
a correlated manner - based on the measured charged- 
particle multiplicity in A+A collisions at one given cms- 
energy y/SNN, and the kinetic and chemical freeze-out 
temperatures (T/,T c h e m) based on the measured identi- 
fied particle pt spectra and multiplicities at the same 
cms-energy. After this, predictions for other ^/snn and 
A are can be computed. 
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The results presented in this paper are quite encour- 
aging for the rather obvious further developments of the 
framework. Following Ref. [l7j . we should next extend 
the study to non-central collisions. This requires a lo- 
calization of the saturation criterion, which in turn calls 
for the spatial dependence of the nPDFs in the colliding 
nuclei. Thanks to the recent developments, these tools 
now exist [43| , so that also the localized saturation study 
can now be brought consistently to NLO. A localized 
saturation will also fix the transverse profiles of the ini- 
tial energy densities in the saturated interior of the sys- 
tem, thus considerably decreasing the profile uncertainty. 
Ultimately, to compute initial conditions for the event- 
by-event hydrodynamics, such a localized study should 
be performed event-by event, using a Monte Carlo sim- 
ulation which accounts for the fluctuations both in the 
number of binary NN collisions and in the minijet mul- 
tiplicity (Et) from one NN collision to another. 

Another very interesting line of further studies is the 
inclusion of dissipation, viscous corrections, to our hydro- 
dynamic description. Dissipation will somewhat increase 
the entropy during the hydrodynamic evolution of the 



system, and the viscous corrections also affect the hadron 
Pt spectra. Promisingly, Fig. Q] indicates that that there 
indeed is room for such an entropy increase, with suit- 
ably chosen values of the parameters /3 and -Ksat ■ It will 
be very interesting to see whether the present framework 
still retains its predictive power also in the presence of 
dissipation. 
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